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We present a detailed description of techniques developed to combine 3D numerical simulations 
and, subsequently, a single black hole close-limit approximation. This method has made it possible 
to compute the first complete waveforms covering the post-orbital dynamics of a binary black hole 
system with the numerical simulation covering the essential non-linear interaction before the close 
limit becomes applicable for the late time dynamics. In order to couple full numerical and perturba- 
tive methods we must address several questions. To determine when close-limit perturbation theory 
is applicable we apply a combination of invariant a priori estimates and a posteriori consistency 
checks of the robustness of our results against exchange of linear and non-linear treatments near 
the interface. Our method begins with a specialized application of standard numerical techniques 
adapted to the presently realistic goal of brief, but accurate simulations. Once the numerically mod- 
eled binary system reaches a regime that can be treated as perturbations of the Kerr spacetime, we 
must approximately relate the numerical coordinates to the perturbative background coordinates. 
We also perform a rotation of a numerically defined tetrad to asymptotically reproduce the tetrad 
required in the perturbative treatment. We can then produce numerical Cauchy data for the close- 
limit evolution in the form of the Weyl scalar ip4 and its time derivative dttpi with both objects 
being first order coordinate and tetrad invariant. The Teukolsky equation in Boyer-Lindquist coor- 
dinates is adopted to further continue the evolution. To illustrate the application of these techniques 
we evolve a single Kerr hole and compute the spurious radiation as a measure of the error of the 
whole procedure. We also briefly discuss the extension of the project to make use of improved full 
numerical evolutions and outline the approach to a full understanding of astrophysical black hole 
binary systems which we can now pursue. 

PACS numbers: 04.25.Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 

Binary black hole mergers are among the most pow- 
erful and efficient sources of gravitational radiation in 
our universe and are thus the primary targets for direct 
experimental detection by the future interferometric ob- 
servatories. Recent astronomical observations of x-ray 
emission sources reinforce the evidence of black holes in 
many galaxies, and astrophysical simulations of globular 
clusters [1, 2] show binary black holes mergers in such 
an abundance to boost the gravitational wave detection 
rate estimation to 1.6 x 10~ 7 yr~ Y Mps~ 3 , which results 
in about one detection event every 2 years for LIGO and 
in one event per day for LIGO II. 

It is thus not surprising that on the theoretical side 
the study of binary black hole mergers has become one 
of the most exciting and challenging topics in astrophys- 
ical relativity. Several theoretical approaches have been 
developed for treating these systems. So far the post- 
Newtonian approximation (PN), has provided a good un- 
derstanding of the early slow adiabatic inspiral, or "far- 
limit", phase of these systems. Similarly, for the final mo- 
ments, when black holes are close enough to each other 
to sit inside a common gravitational well, one can suc- 
cessfully apply the "close limit" approximation (CL) [3], 



which effectively describes the whole system as a pertur- 
bation of a single black hole which rapidly "rings-down" 
to stationarity. Before this last stage, though, when the 
black holes are still close to the innermost stable circular 
orbit (ISCO), the orbital dynamics are expected to yield 
to plunge and coalescence. No approximation method 
can be applied in this highly nonlinear phase and it is 
generally expected that one can only treat the system by 
a full numerical (FN) integration of Einstein's equations. 

Intensive efforts have been underway in the past decade 
to develop numerical codes able to solve Einstein's gen- 
eral relativity equations, by the use of powerful super- 
computers. So far the numerical treatment of black hole 
systems in full three dimensions (3D) has proved very 
difficult and challenging because of the huge computer 
memory requirements, on one hand, and of very severe 
numerical instabilities, on the other, which make the 
codes fail before any useful gravitational wave informa- 
tion can be extracted. In spite of such difficulties, in- 
teresting progress has been made, including, for example 
the work in [4], where a true 3D simulation based on 
the traditional 3 + 1 decomposition of space and time 
has been successfully carried out for the so-called non- 
axisymmetric 'grazing' collisions of two black holes. How- 
ever, because of the limited evolution time achievable be- 
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fore these codes become unstable or otherwise inaccurate 
these simulations must still begin too late in the plunge 
to be practical for direct astrophysical application. In 
most cases treatable so far, the close limit approxima- 
tion theory represents a good alternative model for the 
late time dynamics of these systems. 

Considering the above situation, in Refs [5, 6] we in- 
troduced a new hybrid approach to the binary black hole 
merger problem, called the Lazarus Project, with the mo- 
tivation of providing expectant gravitational wave ob- 
servers with some early estimate of the full merger wave- 
forms within a 'factor two', and to guide future, more 
advanced numerical simulations. The key idea of the 
Lazarus Project is very simple: combine the best of the 
already existing approaches by applying each of these 
methods in sequence and in their best suited regime, 
while focusing the numerical simulations squarely on the 
intermediate phase of the interaction where no available 
perturbative approach is applicable. 

Clearly, the primary task of the combined model is de- 
veloping appropriate interfaces between these three ex- 
isting treatments in such a way that we can also benefit 
of future improvements in any of the above three ap- 
proaches. In an earlier letter [5] we presented the first 
results of our eclectic approach for a model problem, 
the head-on collision of black holes, where we success- 
fully addressed to the problem of combining the close- 
limit approximation describing ringing black holes and 
full three-dimensional numerical relativity. In this well- 
known case, our method proved capable of determining 
radiation waveforms with accuracy comparable to the 
best published 2D numerical results, allowing at the same 
time a more direct physical understanding of the colli- 
sions and indicating clearly when non-linear dynamics 
are important as the final black hole is formed. Pre- 
vious attempts to make a combined use of numerical 
and close-limit evolution [7] have been implemented in the 
case of two axisymmetric black holes formed by collaps- 
ing matter [8], using a 2D numerical code and I = 2 
metric perturbations (d la Zcrilli) of the Schwarzschild 
background and are not generalizable to full 3D simula- 
tions. In Ref.[6] we studied the non- axisymmetric coales- 
cence of equal mass non-spinning binary black holes from 
the innermost stable circular orbit (ISCO) down to the 
final single rotating black hole, and provided the first, 
astrophysically plausible, theoretical predictions for the 
gravitational radiated energy, angular momentum, and 
waveforms to be expected from these systems. 

A sketch of the eclectic approach to the binary black 
hole calculation is outlined in the following steps: (1) 
First provide a description of the early dynamics of the 
system with an approach, such as the post-Newtonian 
method, which is appropriate for slowly moving, well- 
separated black hole. A recent interest within the post- 
Newtonian and gravitational wave research community 
in providing Cauchy data for simulations may soon lead 
to a practical PN-FN interface. (2) Extract critical in- 
formation about the late-time configuration of this sys- 



tem, and translate this information to a corresponding 
solution of the gravitational initial-value problem. (3) 
Apply a full 3D numerical simulation of Einstein's equa- 
tions to generate a numerical spacetime covering the non- 
linear interaction region of the spacetime. The evolution 
should proceed for long enough so that the subsequent 
evolution of the region exterior to the final single rem- 
nant black hole can be well approximated by perturba- 
tive dynamics. (4) At this point we choose a "late-time" 
slice from the numerically generated spacetime, extract 
ip4 = C a fj 1 sn a fh t3 n 1 fh s and d t if>4, to quantify the de- 
viation of the numerical spacetime from a Kerr geome- 
try. Then (5) evolve via the Teukolsky equation, which 
governs the dynamics of Kerr perturbations in the time- 
domain [9] , long enough to drive all significant radiation 
into the radiation zone where it can be interpreted. Mak- 
ing the greatest possible use of perturbation theory in this 
way, not only saves precious three dimensional computa- 
tional resources, concentrating these, for the first time, 
squarely in the intermediate coalescence phase, but also 
provides a new framework to explore and interpret the 
interesting new physics that is expected to take place in 
the transition from nonlinear to linear dynamics. 

The emphasis of this paper is to realize steps (2)- (5) 
above and to describe in detail a general approach to pro- 
viding the FN-CL interface. In Section II we discuss our 
approach to the full numerical simulations which we have 
used to achieve a successful evolution of truly detached 
black holes for the first time. This discussion naturally 
divides into two parts a) our preparation of the initial 
data, by which we greatly improve the simulations effi- 
ciency and b) our numerical evolution method. 

Two important questions arise in implementing the 
transition, step (3), from a numerical approach to a per- 
turbative approach. First, how long must we evolve the 
system numerically before we can obtain a reliable de- 
scription in terms of a single perturbed black hole? We 
use a combination of several independent and comple- 
mentary indicators to establish when perturbation the- 
ory should begin to work. In Section III, we discuss our 
study of two of such indicators: a) The speciality invari- 
ant, S, introduced in Ref. [10], which is exactly equal to 1 
for the Kerr geometry with leading deviations quadratic 
in the gravitational distortions, b) By extracting Cauchy 
data at successive later numerical time slices. When the 
system entered the linear regime, the waveforms evolved 
via the Teukolsky equation should essentially superpose 
to each other. Consequently also a certain level-off of 
the radiated energy should be observed. While a) gives 
a local measure of the physical distortions from the Kerr 
geometry, b) rather depends on the past light cone data. 

The second question is how to identify the single "back- 
ground" black hole which is emerging in the numerical 
spacetime. In order to define deviations from this back- 
ground black hole we must be able to relate it, by an 
explicit diffeomorphism, to the numerical spacetime. We 
need to specify both the spatial coordinates and the time 
slice, which in general may be different from the one 
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used to numerically integrate Einstein equations. This 
geometrical puzzle is discussed in detail in Section IV. 
There is in general no geometrically preferred way to as- 
sociate the numerical and background spacetimes, but 
the first order gauge and tetrad invariance of the per- 
turbative formalism implies that the results should not 
depend strongly on small variations in these choices. 

In Section V, we describe how to compute the Cauchy 
data for the Teukolsky equation, i.e. the Weyl scalar ip4 
and its background time derivative d t ip4, from the nu- 
merical three metric gij and extrinsic curvature K%j, on 
the transition Cauchy hypersurface. The numerical cal- 
culation of the Cauchy data require, first, a nontrivial 
identification of an appropriate numerical 'tetrad', which 
reduces to the (null and complex) tetrad used in the 
perturbative calculation in the small perturbation limit. 
Second, the numerical calculation of dtipi is done 'on 
slice' using Einstein's equations to be consistent with the 
Boyer-Lindquist time of the final Kerr black hole. 

In Section VI we briefly describe the perturbative 
Teukolsky equation and the 2+1 numerical code used to 
numerically to solve it. We then apply all of our tech- 
niques, in Section VII, to evolution of a single Kerr black 
hole with vanishing shift and maximal slicing to test the 
consistency of our method. Our essentially trivial re- 
sult is obtained in a very non-trivial way since our nu- 
merical tetrad is not necessarily aligned with the prin- 
cipal null directions, nor are our numerical coordinates 
the Boyer-Lindquist coordinates used in the perturba- 
tive code. Only after we make the appropriate rotation 
of the tetrad and transform the coordinates to reproduce 
the Boyer-Lindquist ones we see quadratic convergence 
to near vanishing outgoing gravitational radiation. 



II. SUMMARY OF THE FULL NUMERICAL 
TECHNIQUES 

In our full numerical simulations we use many of the 
standard techniques applied in, for example, the Binary 
Black Hole Grand Challenge effort, with adaptations ap- 
propriate the needs of our more specifically defined nu- 
merical simulation problem. Many previous applications 
of numerical relativity to the binary black hole problem 
have been developmental test problems aiming toward 
an ultimate goal of indefinitely long-running 3D numer- 
ical simulations to cover the evolution beginning with 
well separated black holes and evolving through the en- 
tire interaction until further radiation is no longer sig- 
nificant. With regards to gravitational radiation, these 
efforts have been focused on indefinite numerical stability 
and successful radiation waveform extraction by an ob- 
server in the "far away" region of the numerical domain. 
These efforts have often been successful with relatively 
brief black hole evolutions, but have demonstrated the 
serious difficulties in succeeding with the desired long- 
running numerical simulations, and this approach has not 
yet generated radiation studies which approach relevance 




FIG. 1: The eclectic approach: We represent the three phases 
of the binary black hole evolution and the corresponding tech- 
niques adapted to each phase. The full numerical (FN) evolu- 
tion is located to cover the truly nonlinear dynamical interac- 
tion. The domain of perturbative evolution (CL) follows the 
FN domain allowing indefinite evolution. Waveforms are ex- 
tracted at the dotted world line depicted on the right. Though 
such observers are located in the CL part of the spacetime 
they will experience all radiation arriving from the strong field 
dynamical FN region. In the far limit regime we envision to 
use the post-Newtonian (PN) approximation 



to astrophysical problems. 

We will ask less of our numerical simulations. Our de- 
mand is for a highly accurate determination of the most 
significantly non- linear part of the binary interaction. We 
will try to make use of codes that may only run stably 
for a relatively brief period, but which can provide an 
accurate representation of the part of the spacetime we 
are most interested in. This point of view allows us, 
for example, to avoid the difficult problem of imposing 
physically accurate outer boundary conditions, by only 
considering the part of the spacetime causally separated 
from the boundary. We find that this can be done much 
more efficiently in specialized coordinates, described in 
the first section below. Similarly, we have not yet needed 
more stable formulations of Einstein's equations, or diffi- 
cult sophisticated techniques such as black hole excision. 
Our straightforward numerical approach to evolution is 
described in the second section. 



A. Preparing the initial data 

Ultimately we wish to derive initial data based on in- 
formation from an approximation procedure, such as the 
post-Newtonian method which is applicable in the limit 
of slow-moving/far-apart black holes. As no such inter- 
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face is presently available we use, in our present work, 
initial data from an alternative source, commonly ap- 
plied in numerical relativity, the "puncture" formalism 
with conformally flat three-metric and purely longitudi- 
nal extrinsic curvature on a maximal slice. This assumes 
a three-sheeted topology instead of an inversion symme- 
try across the throats [11] allowing for a solution of the 
elliptic Hamiltonian constraint equation without having 
to impose interior boundary conditions [12]. 

Within this family it is possible to identify data 
roughly corresponding to quasi-circular orbits using the 
effective potential method as in Ref. [11]. The binding 
energy of the system is computed as a function of the 
proper separation of the holes keeping everything else 
constant. A minimum in the binding energy is then in- 
terpreted as giving an stable quasi-circular orbit. Within 
this approach an ISCO is determined by varying the or- 
bital angular momentum of the system until this min- 
imum becomes an inflection point. For less separated 
configurations, a stable quasi-circular orbit is no longer 
possible. We use these ISCO data, determined in [13], 
for non-spinning equal-mass black holes as a particularly 
reasonable starting point for approaching astrophysical 
systems [6]. 

Having selected the physical initial data we then pre- 
pare it for numerical evolution. When Smarr and York 
[14] spelled out the problem of 3+1 numerical relativity 
in the 1970's, they specifically sought out methods which 
would be invariant to gauge transformations in the initial 
data. In the pursuit of long-running all-purpose numeri- 
cal relativity tools, this viewpoint has been traditionally 
preserved, and little attention has been given to the ques- 
tion of choosing appropriate coordinates for the initial 
data. It is clear though, that whenever differential equa- 
tions are to be solved numerically, some choices of vari- 
ables (coordinates) will be more practical than others. In 
a wave-propagation problem, for instance, the simulation 
will be much more efficient if a wave is evenly resolved as 
it moves across the numerical domain, or similarly, if the 
coordinate characteristic speeds were constant in space 
and time. 

For numerical relativity simulations in practice we are 
often very far from this ideal. In typical coordinates, 
such as isotropic coordinates for our (initially) confor- 
mally flat spaces, the waves are strongly red-shifted as 
they move away from the strong-field region. Since we re- 
quire both a physically large computational domain and 
also high resolution in the strong field region, use of the 
standard coordinates leads to a great waste of numerical 
effort on over-resolving an outgoing radiation wave which 
was originally generated with much poorer resolution. In 
this way, relatively little is gained by, expanding the com- 
putational domain with additional numerical grid-points. 
We find that we can make great improvements in numer- 
ical efficiency with a relatively simple ad hoc coordinate 
transformation on the initial data which we call 'fish-eye' 
coordinates. A typical such transformation is a radial 
rescaling, r lso = R num cosh ((R num / R ) n ) with typical 



values Rq = 7.7 and n — 2. This allows us to maintain 
a central resolution of up to M/24 with outer bound- 
aries near r new = 37 M using only 256 x 512 2 grid-points, 
moving the outer boundary much farther away without 
loss of physical resolution in the strong field region. This 
problem is illustrated in Figure 2, which shows data from 
numerical simulations in two alternative coordinate sys- 
tems after 10M of evolution from an initial ISCO config- 
uration. The outgoing radiation wave is noticeable in the 
real part of the gauge-independent S invariant discussed 
in Sec. Ill A. These curves represent the same physical 
spacetime as seen from alternative numerical coordinate 
systems. In this figure the strong field dynamics are most 
important in the left side on the figure up to about the 
value of the numerical coordinate (along the z-axis) of 
Znum — 6. Up to that point the two coordinate sys- 
tems are nearly identical. As we add grid-points on the 
right side of the figure beyond this strong field region, 
we are frustrated, in the isotropic coordinate case by the 
red-shift effect and only modest additional part of the 
outgoing wave, about half a wave cycle, is added to the 
grid when we roughly triple the grid dimension. In the 
case of our fish-eye coordinate the wave is evidently more 
evenly resolved across the grid, and we cover the domain 
of the isotropic coordinate system with only about a 60% 
increase in the grid dimension. As shown, the fish-eye co- 
ordinate system has a much more distant physical outer 
boundary than that of the isotropic case while still hav- 
ing only about half of the grid-points in 3D. Note that 
we would gain no further advantage by attempting to 
compactify spatial infinity as in for example [15] since 
resolution must nevertheless still fail to resolve waves at 
a finite radius in such a scheme. 

Foreseeing longer term full numerical evolutions we 
have also implemented other rc-coordinatizations of the 
initial data that have a fairly constant high resolution in 
the center of the grid (where the grid stretching is more 
severe) and a lower resolution near the boundaries, but 
still fairly constant to allow the application of the usual 
radiative boundary conditions (adapted to the different 
characteristic speed) . One of such functions is 

r iso = R num (l + b (tenh ^ R ^yRo) + 35 ^J 

+ tanh(^-0.35^ (2.1) 

with b, d, Rq adjustable parameters that determine the 
ratio of central to boundary resolutions, the width and 
location of the effective resolution transition region re- 
spectively 



B. Numerical evolution 

Our numerical evolution must be consistent with our 
need for highly accurate relatively brief simulations. 
Consequently, in our work so far, we have used the stan- 
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FIG. 2: The benefit of our fish-eye coordinates compared 
against the typical isotropic coordinates. The S invariant, 
plotted here, gives an indication of the radiation moving out 
from an initial ISCO system after 10M of numerical evolution. 
In the strong field region up to z = 6 the two coordinate sys- 
tems are very similar. Moving outside that region though the 
Fish-eye coordinate cover a significantly larger region of the 
physical spacetime with fewer grid-points. The extra grid- 
points in isotropic coordinates are wasted by over-resolving 
the outer part of the radiation. In Fish-eye coordinates the 
wave is resolved more evenly. 

dard ADM (Arnowitt-Deser-Misner) formulation of Ein- 
stein equations [16] as adapted by Smarr and York [14]. 
Our evolution equations are thus simply: 

Bodab = -2aK ab (2.2) 
d K ab = -V a V 6 a 

+a(R ab - 2K ac K c b + K ab K) (2.3) 

where do — dt — £p. Here and below Latin indices run 
from 1 to 3. 

Though a newer conformal formulation of Einstein's 
equations has been found to be more stable in various 
numerical simulations [17], here we focus on the accu- 
racy of the solutions rather than long term stability. Our 
observation is that the standard ADM equations seem to 
give more accurate results for binary black hole simula- 
tions in our gauge while the simulation is stable. 

If it is possible to have a slicing which is consistent 
with that of our perturbation theory, then we can avoid 
a rather large technical problem of producing data on a 
slice inconsistent with the background. Consistent with 
our choice of Boyer-Lindquist coordinates in our pertur- 
bation treatment of the background black hole, we have 
chosen maximal slicing to define the lapse a, 

K = 0, Aa = a K ab K ab . (2.4) 

This implies an elliptic equation for a which we typically 
have solved every 5 time steps using Dirichlet bound- 



ary conditions. For simplicity we set the shift /3 ! = 0, 
which is an adequate condition for relatively brief runs. 
The numerical evolution is performed using an iterative 
Crank-Nicholson method of third order which is second 
order convergent. In our simulations we have used a res- 
olutions up to dx = M/24 with dt — 0.25 dx. Because we 
have moved the outer boundary to a point causally sepa- 
rated from the region we are interested in it is acceptable 
simply to impose static boundary conditions. 

In evaluating the results of our numerical simulations 
we make frequent use of two indicators: The degree of 
satisfaction of the ADM constraint equations gives a mea- 
sure of the numerical error produced by the evolution 

V a (i^ ab - g ab K) = (2.5) 
R - 2K ab K ab + K 2 = 0. (2.6) 

These quantities provide an important indication of when 
numerical inaccuracies (and eventually instabilities) have 
become significant in our simulations. Even if Einstein's 
equations could be solved perfectly, any simulation with 
a finite boundary is subject to an additional type of er- 
ror arising from inappropriate boundary conditions. A 
geometrically correct solution may have physically un- 
reasonable disturbances propagating in from the bound- 
ary. We have found the speciality invariant S [10] to be 
a sensitive indicator such boundary waves, which do not 
violate the constraints. 



III. DETERMINING THE LINEAR REGIME 

Black hole perturbation theory has recently generated 
much interest as a model for the late stages of a binary 
black hole collision spacetime [3] . When two black holes 
are close enough to each other one can simply treat the 
problem, in the 'close limit' approximation, as a single 
distorted black hole that 'rings down' into its final equi- 
librium state. So after some nonlinear numerical evolu- 
tion of the full Einstein's equations for a system of two 
initially well-detached black holes, there should always 
be a transition time, T, after which the system sim- 
ply behaves linearly i.e. satisfy the linear perturbation 
equations around the final Kerr black hole. Finding the 
linearization time T is thus the first nontrivial question 
which arises in the context of our 'eclectic' approach. 
In other words, we need one or more working criteria 
for when we can expect perturbation theory to be accu- 
rately effective based only on numerical data. As we shall 
see below, we apply at least two independent criteria for 
estimating the onset of linear dynamics, the speciality 
invariant prediction based only on the Cauchy data and 
another estimate based on the stability of the radiation 
waveform phase. 
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A. The speciality invariant test 

Motivated by this purpose in Ref. [10] we introduced 
an invariant quantity, 



Re(S), along z-axis 



5 = 27J 2 /! 3 , 



(3.1) 



where X and J are the two complex curvature invariants 
X and J , which are essentially the square and cube the 
self-dual part, C a()l s = C a/3l s + {i/2)e abmn C m ^ d , of the 
Weyl tensor: 



X — CafjysC ^ and CI — C a p^$C^ ^ V C^ 



(3.2) 



Both these scalars can be expressed in terms of the Weyl 
components, for an arbitrary tetrad choice: 

X = 3ip 2 2 - 4V>i^3 + -04^0, 

J = + ipQtp A rp 2 + 2ipiip 3 ip 2 ~ - V^f-3-3) 

The geometrical significance of 5 is that it measures 
the deviations from algebraic speciality (in the Petrov 
classification of the Weyl tensor). 

For the unperturbed algebraically special (Petrov type 
D) Kerr background 5 = 1. However, for interesting 
spacetimes involving nontrivial dynamics, like distorted 
black holes, which are in general not algebraically spe- 
cial (Petrov type I), we expect 5 = 1 + A5, and the size 
of the deviation AS ^ can be used as a guide to pre- 
dict the applicability of black hole perturbation theory. 
In particular we adopt the criterion that, when 5 differs 
from its background value of unity by less than "a factor 
of two" outside the (background) horizon, a perturbative 
treatment may be expected to provide a reasonable de- 
scription of the radiative dynamics. A larger deviation 
from algebraic speciality implies significant "second or- 
der" perturbations. In fact, for perturbations on a back- 
ground Kerr spacetime, with an arbitrary tetrad pertur- 
bation, one can easily deduce 



5 = 1 - 3e 



(3.4) 



where ipo , tp4 and ip2 are the usual Newman-Penrose com- 
plex Weyl scalars. The lowest order term in the devia- 
tion is second order in the perturbation parameter e, and 
should tend to vanish if first order perturbation theory 
is appropriate. Note that the superscripts (0) and (1) 
stand respectively for background and first order pieces 
of a perturbed quantity, where e is a perturbation pa- 
rameter. 

In Fig. 3 we display the speciality invariant along the 
z-axis, perpendicular to the orbital plane of two black 
holes starting the evolution from the ISCO determination 
used in[6]. Its value oscillates around one (the Kerr back- 
ground value). After some evolution T w 11M, the am- 
plitude of the oscillation decreases to a deviation below 
50% outside the horizon (located at around Zjjtjm ~ 2.5 
in the numerical coordinates, and perturbation theory 
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FIG. 3: The 'speciality' invariant for binary black holes evolv- 
ing from the 'ISCO' showing damped oscillations around 
unity, its Kerr value. The location of the horizon in these 
coordinates is roughly 2.5. Its behavior at larger radius sug- 
gests radiation is beginning to leave the system. 



can reliably take over the remaining of the evolution. 
Because the gravitational field has two degrees of free- 
dom is it clear that the 5-invariant alone is insufficient 
to provide a complete description of black hole perturba- 
tions, and can be complemented with its time variation 
5. Consequently, we have been looking at the turning 
points where 5 = and the amplitude of the distortion 
reaches a maximum. 

As noted in Section III, the 5 is also very useful out- 
side the perturbative context. Its usefulness is derived 
from the fact that it is a gauge invariant quantity which, 
unlike / and J is not dominated by strong "peeling prop- 
erty" fall-off behavior, which tends to indicate spacetime 
dynamics. Because the Weyl tensor, C a /3~fS, carries infor- 
mation about the gravitational fields in the spacetime, 5 
turns out to be an interesting indicator of radiation of 
the spacetime and tests, for instance, how much radiation 
is produced by the imposition of approximate boundary 
conditions. We have found that the 5— invariant is simple 
to calculate and can be applied directly to full 3D numer- 
ical evolutions to provide a gauge invariant indication of 
the dynamics. 



B. Waveform locking and energy plateau 

The phase and the amplitude of the radiation, or equiv- 
alently the locking of the waveforms and the correspond- 
ing energy plateau, also provide an indicator of linear dy- 
namics. Starting with detached black holes, we expect an 
initial period of weak bremsstrahlung radiation followed 
by the appearance of quasi-normal ringing. On the other 
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ISCO, total energy radiated RefyJ , r*=30M, 9=0 




FIG. 4: Energy radiated from two black holes from ISCO ^ \ j 

configuration for different transition times showing a plateau ~® ^30 35 40 45 50 55 

when reaching the 'linear' regime. t/M 



hand, switching to perturbative evolution immediately 
leads to premature ringing. Hence if we cut short the 
numerical simulation and apply linear theory too early, 
we observe quasi-normal ringing too early and calculate 
a waveform which is out of phase with the desired result. 
Comparing waveforms derived from differing durations of 
numerical simulation then we tend to see a phase shift 
in the onset of the ringing when we have not yet allowed 
enough numerical simulation. In practice, we thus follow 
the behavior of the waveforms through the evolution by 
extracting the Cauchy data at successive later numerical 
time slices. When the system enters in the linear regime, 
the waveforms evolved via the perturbative Teukolsky 
equation should essentially superpose to each other, as 
changing this transition time amounts to an equivalent 
exchange of linear and non-linear evolution for the inter- 
vening region of spacetime. Consequently also a certain 
level-off of the radiated energy should be observed (See 
Fig. 4). 

As we show in Fig. 5, extracting waveforms every 1M 
of non-linear numerical evolution allows us to study the 
transition to linear dynamics, and to perform important 
consistency tests on our results. If we have made a good 
definition of the perturbative background, as described 
in Section IV, then we can expect our radiation wave- 
form results to be independent of the transition time, T, 
once the linear regime is reached and for as long as the 
numerical simulation continues to be accurate. 

A closer look at Fig. 5 gives us an idea of how the lin- 
earization happens. Curves of T = 10&11M of evolution 
are close to the correct waveform for this orbital case 
starting at a proper separation L/M — 4.9. If we apply 
right away the close limit approximation we get the curve 
labeled by T — 0M which starts ringing prematurely. Af- 
ter 2M of full numerical evolution we obtain good agree- 
ment with the correct waveform up to t/M w 33. When 
perturbation theory takes over after AM of full numerical 
evolution the agreement is very good up to t/M « 38. 



FIG. 5: Detail of the progressive waveform locking process for 
black holes at the location of the ISCO. 

Near t/M = 45M we seen we need 8M of nonlinear evo- 
lution while near t/M = BOM 10M of full numerical 
evolution are needed and at longer times the agreements 
begins to be fine for the whole relevant waveform. This 
process shows how the full nonlinear dynamics shifts to 
a central region covered by the common potential barrier 
allowing to describe linearly the evolution of the outer 
part. 

C. Common horizon 

An intuitive picture to visualize the applicability of 
the close limit approximation would be the appearance 
of a common event horizon that encompasses the binary 
system. Under these conditions the spacetime exterior 
to the horizon (the relevant one for computing gravita- 
tional radiation reaching infinity) can be treated as per- 
turbations of a Kerr hole. In practice event horizons are 
difficult to compute in numerical relativity because they 
are a global feature of the spacetime and we would need 
to first evolve the binary system for a long time and then 
extract a posteriori the information to locate the event 
horizon. An easier quantity to compute is the appar- 
ent horizon that can be defined locally as the outermost 
marginally trapped surface of the spacetime where a con- 
gruence of null rays directed outwards have vanishing ex- 
pansion [18]. A common apparent horizon lies inside and 
in a binary system appears later than a common event 
horizon; and typically much later than when the sys- 
tem can be effectively described by linear perturbations. 
The linearization time refers to when the close limit ap- 
proximation can be applied and this intuitively occurs 
when a common potential barrier covers the binary sys- 
tem. In black hole perturbation theory, a potential is 
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present somewhat outside the horizon of the black hole 
which tends to prevent radiation from escaping this re- 
gion. This is the main reason why the close limit is such 
a good approximation even beyond original expectations 
[19]. 



IV. CONSTRUCTING THE KERR 
BACKGROUND 



(t, r, 8, <fr) reads, 



ds 2 



2Mr\ , 2 /XT, , 
" dt 2 + f - ) dr 2 



+zZd6 2 + sin 2 6%d<j> 2 - sin 2 6dtd<t>, (4.1) 



where A 



2Mr 



,S = r 2 + a 2 cos 2 9 and 



(r 2 + a 2 )S + 2Mra 2 sin 2 6>, M is the mass of the black 
hole, a its angular momentum per unit mass. 



The slice 



Einstein's theory of gravity in principle demands the 
equivalence of all coordinate representations of gravita- 
tional dynamics. However, in practice one always needs 
to choose a convenient gauge to accurately carry over the 
full numerical evolution to the point where the two black 
hole system effectively behaves like a single perturbed 
black hole. Having determined that a late time numerical 
spacetime geometry is close to the Kerr spacetime does 
not give us any information about the coordinate system 
in which this is written. In order to be able to con- 
tinue the numerical evolution with the Tcukolsky equa- 
tion (6.2), we thus need to reconstruct a Kerr background 
in a recognizable form, for instance in Boyer-Lindquist 
coordinates. Because there is in general no unique pro- 
cedure to reconstruct such a Kerr background, we shall 
require that this should be close enough to the given nu- 
merical spacetime. In other words, we will require that 
the two spacetimes agree to the first order in e. Since the 
physics of our problem will then be described by quan- 
tities, like "04 j which are first order gauge (and tetrad) 
invariant, the physical results we compute will be inde- 
pendent (to first order) of the identification of the back- 
ground coordinates we describe below. To have complete 
theoretical control of the perturbation theory, it is de- 
sirable to have to a complete family of initial data sets 
which reduces to the background geometry in the limit 
e — > 0. While this requirement is not strictly required in 
a practical perturbative application [7], we would like to 
stay as close as possible to this arrangement for its bene- 
fit in evaluating our results. In our case the perturbation 
parameter e can be regarded as a decreasing function of 
the transition time T. In practice, we will not be able 
to achieve an exact Kerr black hole in the T — > oo limit, 
but we will aim for the practical goal that the remaining 
perturbations are small compared to the radiation we are 
interested in, a condition which we test in Section VII. 

We initially suppose that the background Kerr black 
hole is given by the parameters M and a of the initial 
data. With a first estimate of the total radiated energy 
and angular momentum these parameters can be iterated 
to approach the final values for the stationary Kerr black 
hole. 

The Kerr metric in Boyer-Lindquist coordinates 



We recall that a Boyer-Lindquist slice of the Kerr met- 
ric has K — 0. The full numerical coordinate condition 
of Maximal slicing, Eq. (2.4), is solved for the lapse a 
with an exterior boundary condition set to reproduce 
the value of the Boyer-Lindquist lapse there, but to van- 
ish at the location of the individual black hole "punc- 
tures". The resulting lapse from the evolution of two 
holes from 'ISCO' is shown in Fig. 6. The lapse re- 
sembles the Boyer-Lindquist lapse initially and further 
evolution bring them closer. Thus the maximal lapse 
with our boundary condition approaches the background 
lapse quite closely. Where there are differences, near the 
horizon, our lapse tends to produce a coordinate system 
in which the coordinate observers drift slowly into the 
black hole. Considering our coordinate trajectories from 
the frame of the background black hole, one can conclude 
that since the trajectories and lapse are similar away from 
the horizon, and the lapse is a bit different near the hori- 
zon, our slicing will be close to the background slicing, 
but slightly distorted toward the future near the horizon. 
In Section VI we try to quantify the significance of this 
distortion with a numerical study of the Kerr spacetime 
in these coordinates. 

Other lapse possibilities can be considered which also 
produce a slicing similar to that of the Boyer-Lindquist 
background, algebraic slicings, for instance [20] l (l+log)\ 
and a re-parametrization of the maximal slicing by an 
f(a) such that the numerical lapse resembles even closer 
the Boyer-Lindquist one. We performed such tests and 
check that whenever the deviations from the Boyer- 
Lindquist lapse are close enough the results for radiated 
waveforms and energies do not change notably, in agree- 
ment with the first order gauge invariance of ip4 and 9^4- 

B. The spatial coordinates 

The general idea here is to numerically compute phys- 
ical quantities or geometrical invariants and relate them 
to their analytic expressions in the pcrturbativcly pre- 
ferred coordinate system. Curvature invariant methods 
have the distinct advantage that they can be applied to 
evolutions using numerically generated coordinates which 
are not understood analytically. On the other hand, the 
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Maximal lapse, along z-axis 




r*/M 



FIG. 6: The maximal lapse used for black holes evolving 
from 'ISCO' compared to the analytic Kerr lapse in Boyer- 
Lindquist coordinates for different evolution times. 



The second important coordinate effect, which be- 
comes relevant when the total angular momentum is sig- 
nificant, is the result of frame dragging caused essentially 
by the difference between our vanishing shift, and the 
non- vanishing Boyer-Lindquist shift. This effect drags 
the coordinates in the ip direction and has the effect 
of producing an off-diagonal distortion in the numerical 
metric. We can undo the frame dragging by attempting 
to restore the diagonal form of the Boyer-Lindquist three 
metric. 

We seek to set frame-dragging gauge freedom by sup- 
plementing the Cartesian definition of ip = arctan[y/x] 
with a correction that makes the metric component most 
strongly affected by frame dragging, g r 0, vanish, 

(j) = arctan[y/a,] + / \g rip / g vv )dr, (4.4) 

where the 'hat' stands for the full numerically evolved 
metric. To see this, consider g a b with no shift, and trans- 
form to g ab with ^-shift by p — > 4> — p + tp Eset(t,r,9) 
(note that g^ = j w ). Then, since the Kerr three-metric 
is diagonal, 



so that 



g r <t> = = g rip - d r ip ff BCt g vv (4.5) 

OrPottsct = -T^- (4.6) 

9w 



Similarly, since the numerical metric has zero shift, we 
find 



values of curvature invariants in the perturbed space- 
time may be sensitive to perturbative distortions, making 
them less useful for identifying a background spacetime. 
In light of these effects, we pursue a combined approach, 
utilizing both gauge and geometrical information where 
each seems most appropriate. In the outer regions of our 
spatial slices we expect the gauge to be close to the quasi- 
isotropic gauge for Kerr data. Moving in from this to the 
interior region we expect, most importantly, two gauge 
effects. First, our slicing has the tendency (without a 
shift) to cause the coordinates to fall inward with evolu- 
tion. We counteract this with a rescaling of the radius gtj, -^Kerr 9<P<P 
^Kcrr = ^KcrrM, making use of the X invariant which de- tSPofiset — * — ~ 
pends most significantly on the radial coordinate in the 

background slice, I = 3M 2 /(r - iacos6>) 6 . We use this Since N$_ eII is constant in t, 



5*0 = 9ttp — dt (^offset 9<ptp — —dtPoSset (jtptp 

so that 

JV.t OJJ , 

-iVl„. 



relation and information about the numerical value of X 
to define the rescaled radius. To do this we need to pro- 
duce one value of "Z" for each constant r sphere in the 
numerical slice. The maximum value of I dp tends 
to lie on the equatorial symmetry plane of our binary 
black hole problem, where the 5-invariant also indicates 
relatively weaker distortions. This makes 



(^offset = 0, 



^offset 



-tN, 



Kerr ' 



(4.7) 



(4.8) 



(4.9) 



(4.10) 



1 

2^ 



271 



<1> = 

r K err = {/3A// < 1 > 







l{r,0 = n/2,p)dp (4.2) 



(4.3) 



a practical definition which counteracts the coordinate 
infall. 

Unlike for r there are no obvious dynamical effects on 
the 9 coordinate and it has been sufficient to adopt the 
numerical value cos 9 = zj \j (x 2 + y 2 + z 2 ). We success- 
fully applied these re-mapping of coordinates already in 
the head-on collision case [5]. 



Equation (4.10) allows us to test how close our derived 
(from the block diagonal metric condition) shift correc- 
tion is to the Boyer-Lindquist shift. The results of this 
comparison are displayed in Fig. 7. For two black holes 
evolving from the ISCO, the shift correction correctly re- 
produces the frame dragging effect outside the potential 
barrier of the system and evolution bring the shift closer 
to that of a single rotating Kerr hole. 

We note that some means of fixing this frame-dragging 
degree of gauge freedom, as we have done here, is essen- 
tial also if one wishes to speak meaningfully of the num- 
ber of orbits the system has undergone in the strong field 
region during numerical simulations. 
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FIG. 7: The effective shift correction for black holes evolving 
from 'ISCO' compared to the analytic Kerr shift in Boyer- 
Lindquist coordinates for successive evolution times. 



As already pointed out there is no unique way of choos- 
ing the coordinate transformations in order to bring them 
closer to that of the Kerr background. Our philosophy 
in this section has been to consider the simplest of these 
transformations that approaches the Boyer-Lindquist co- 
ordinates with enough accuracy for the binary black hole 
numerical simulations we are interested in. Obviously, 
other possibilities that would improve the accuracy of 
the procedure can be incorporated as needed. We also 
note that the optimal choice of coordinate transforma- 
tion needed here may depend on the shift condition used 
in evolution and the coordinates used for the initial data. 
The use of a shift condition, such as minimal distor- 
tion (with the appropriate boundary conditions), which 
is naturally adapted to the stationarity Killing vector of 
the background Kerr spacetime [14] may, for example, 
eliminate frame dragging and thus reduce the need for a 
transformation such as (4.4). 



CONSTRUCTING THE CAUCHY DATA 



^4 = [Rijkl + 2K l[k K l]:j ] n l m 3 n k m l 



K 



j[k,l] 



' TP j[k K l\P 



+4 [Rji - K jp Kf + KKji] n [0 m j] n [0 m l \ (5.1) 
and its time derivative 



-N% (^4) 



d R 



i /L I 



n l m^n k m l 



d K m] +d T p jlk K l]p + T p . lk d K l]p 



d Rji - 2K^d Kj) p 



2NK&KPK? 



+Kjid K + Kd Kji n [0 fh j] n [0 fh l] 

+2{ip A (k& - m^N* + 4> 3 (m8 - miA)JV*}. (5.2) 

where the last term extends the expression in the refer- 
ences, having been added to take into account the vari- 
ation of the tetrad terms do [nmnm] . Here A = n^d^ , 
$ — fh^dfj,, and 

ip3 = [R t3 ki + 2K l[k K l]j ] l l n j m k n l 



K 



(l^vP ] fh k n l - n^ Q mP ] l k n l ) (5.3) 
+2 [Rjt - KjpKf + KKji] (^V'JmV - l^n°m l ), 

where the background (null and complex) tetrad, 
(l^, n^ 1 , m M , to m ) is specified in the subsection below. 

The derivatives involved in the above expressions can 
be computed in terms of the initial data on the Cauchy 
hypersurface as in Eq. (A19) in the Appendix. 

With the tetrad specified, the foregoing formulae are 
coordinate independent. Therefore the only adjustment 
needed to specify initial data for the evolution equations 
we will be to insert the appropriate background quanti- 
ties in the above equations. In particular, taking and 
N % respec tively a s the zeroth order Kerr lapse and shift, 
A (0) = v/AE/ft and N l {0) = [0, 0, -2aMr/Q], allows us 
to compute dtip4 directly with respect to the background 
Boyer-Lindquist time, thus avoiding additional perturba- 
tions introduced if one computes the numerical derivative 
by finite differences of '04 on two successive slices. 



Given the numerical metric gij and the extrinsic cur- 
vature Kij derived as in Section II on a Cauchy hyper- 
surface, and the coordinates of the background metric 
determined in Section IV, we proceed to compute the 
Weyl scalar ip4 and its background time derivative 9* ^4, 
the Cauchy data we will need to continue the evolution 
via the Teukolsky equation. As was discussed in refer- 
ences [21, 22, 23, 24], one can make the following 3+1 
decomposition, using the basis 9° — dt, 6 l — dx % + N z dt, 



A. The tetrad 

A null and complex 'exact' tetrad (i.e. orthonormal 
in the numerical spacetime) must be chosen such that 
it reduces, in the linear regime to the choice made in 
our perturbation treatment of the final Kerr hole, the 
Kinnersley tetrad [9] . In Boyer-Lindquist coordinates the 
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background tetrad vectors are 



'Ki 



TV- 



Kin 



l[(r 2 +a 2 ),A,0,a], 
^ [(r 2 +a 2 ),-A,0,a] , 



-\/2(r + za cos 0) 



zasin#, 0, 1, 



sin6 



(5.4a) 
(5.4b) 
(5.4c) 



The Kinnersley tetrad is particularly well suited for per- 
turbation studies because it has the property that and 
are chosen to lie along the (background) principal 
null directions of the Weyl tensor (PND) in such a way 
that one can derive decoupled perturbation equations. In 
terms of the 3+1 basis of Eqs. (5.1)-(5.3), we have 



-^(O)^Kini ^Kin + ^(O)^Kin 
^(0)<i„,<i„+N(0)<: 

A (0) r< in , m l Kin + N l {0) m^ in 



(5.5a) 
(5.5b) 
(5.5c) 



To numerically determine an 'exact' tetrad we could in 
principle search for possible candidates of the two of the 
PND of the Weyl tensor. One could of course, try to pick 
up some null directions in our numerical spacetime g™ m 
which we know are close to the PND in Kerr whenever 
g™" m is a perturbation of Kerr. However, this turns out 
to be a bad choice because the PND do not behave ana- 
lytically under analytic perturbations of Kerr. The rea- 
son is that the principal null directions of Kerr are double 
principal null directions of the Weyl tensor, which in gen- 
eral will split under the perturbation. It turns out that 
the splitting of eigenvectors of an endomorphism under a 
perturbation of order e behaves in general as some frac- 
tional power of e (hence non-smoothly). So, the principal 
null directions will be too strongly perturbed. 

An alternative and more effective procedure to define 
an exact tetrad that has the required property in the lin- 
ear regime is the following, (a) We assume the following 
3 + 1 decomposition of the tetrad 



V 2 

m^ = -L(r+^), 



(5.6a) 
(5.6b) 
(5.6c) 



where is the normalized time-like unit normal to the 
hypersurface and = [0,v%],9» = [0,v§],ipi* = [0,<] 
are orthonormal vectors pointing along the numerically 
defined coordinate directions, (b) We thus identify the 
set of null rotations to bring (5.6) to the form (5.5), in 
order to make it consistent with the tetrad assumed in 
the perturbative calculation. 

Step a is straightforward. Begin with real vectors 
aligned with the numerical space's <p and radial direc- 



tions, which in Cartesian coordinates read 
< = [-y,x,0], 



"2 = [x,y,z] , 

& A / \ 1/2 ad b 

v 3 = dct(g) 1 g e dbc v 1 



b„.c 



(5.7) 



We then redefine these, to achieve ortho-normalization. 
It is important to begin the ortho-normalization proce- 
dure with the azimuthal direction vector v® which is not 
affected by the frame-dragging effect discussed in Section 
IV B. At each step, a Gram-Schmidt procedure is then 
used to ensure that the triad remains orthonormal, so 
that 



«2 



/u>u 



where w tJ = (vfv^gab)- In the case of Kerr one finds, in 
Boyer-Lindquist coordinates, 



1,0,0, 



2aMr 



0M 



[o,<] = 




0,0,-^,0 
0,0,0, 



sin 9 



(5.8a) 
(5.8b) 
(5.8c) 
(5.8d) 



normalized such that —u^u^ = r^r* 1 = 9^9^ = (p^f^ = 1 
so that the inverse metric can be expressed as g^ v = 
2(n>m^ -l^n^) 

For step b identify a combination of null rotations of 
type I and II parameterized by A, and a type III (boost) 
null rotation parameterized by Fa and Fg which bring 
the orthonormal tetrad (5.6) to the form (5.5) for the 
unperturbed case. The transformation 



Fa 
2 

-iA(fh' 1 



(VA 2 + 1 + 1) f" + {y/A 2 + 1-1) fi^ 



')]■ 



(5.9a) 



mr = 



{VA 2 + 1 - 1) P + {yjA 2 + 1 + 1) ii» 

(5.9b) 

(sjA 2 + 1 + 1) - (VA 2 + 1-1) 



-iA{m» - jft")] , 
F B 



+iA(l» + n») 



(5.9c) 



achieves this with A — asm9^ A/fi, Fa — ^/2S/A and 
Fb = (r + ia cos 9) , thereby producing a tetrad con- 
sistent with the tetrad assumed in the perturbative cal- 
culation. 
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In practice we perform the tetrad transformation indi- 
rectly , implementing its effect on the set of Weyl scalars 
(tpo, ... ,^4) as described in the Appendix, Eq. (A24). 



VI. THE TEUKOLSKY EQUATION 



the Lax-Wendroff algorithm, using the standard tortoise 
coordinate r* , 



r = r 



In 



r+ 



2M 



r±=M± \/M 2 - a 2 



In 



2M 



(6.3) 



Perturbations of a rotating Kerr black hole are de- 
scribed by the well known Teukolsky equation [9] , which 
is derived from the Newman-Penrose formalism. The 
Weyl scalar t/> 4 that represents outgoing gravitational ra- 
diation satisfies a decoupled wave equation 

\ (A + Ap + p + 3 7 - 7) (D + 4e - p) 



(S + 3a + j3 + 4n-T)(5 + 4:p-T) 



0. 



(6.1) 



In this generic form the Teukolsky equation is manifestly 
independent of the choice of coordinate system used to 
describe the Kerr background and its perturbations. In 
the foregoing equation the usual notation for spin coeffi- 
cients a, (3, ... was used and A = n^d^, S 
D = l^dn represent directional derivatives. 

For the applications in this paper we consider Boyer- 
Lindquist coordinates (t, r, 9, <p) and the Kinncrslcy 
tetrad. The Teukolsky equation then reads 



m^d^, and 



(r 2 +a 2 Y 



A 



2 • 2 , 
a sin 1 



9 2 V> 4Mar d 2 ip 
W + A dtdcf) 



+ 



1 



a 

A ~ sin 2 6» 
1 d 



d 2 jj 



^0 09 \" m "% ] 



dr \A dr 
M(r 2 - a 2 ) 



which naturally leads to excision of the black hole in- 
terior and constant characteristic wave speed. We im- 
pose static boundary conditions on the internal bound- 
ary (event horizon of the Kerr background) and radiative 
boundary conditions on the exterior boundary. Frame 
dragging effects are taken care of by the background 
Boyer-Lindquist shift. Thus, this formulation has all the 
ingredients to allow for an indefinite stable evolution. In 
practice it provides an accurate evolution for the few hun- 
dreds of M of relevant signal generated in the final stages 
of black hole merger. Since the Kerr background has the 
axial killing vector we can Fourier decompose ip4 into 
e im4>Ks mo des. In particular, for numerical convenience, 
we use the 'Kerr-Schild' 4>ks 



4>KS 



+ 



a 



■In 



r+ — r 
a 



r + — r_ 



In 



±-1 



(6.4) 



A 



r — ia cos 9 



~dt 



+4 



a(r - M) 



l cos( 



sin 9 



^ + (4cot 2 (9 + 2)?A = 0, (6.2; 



where ip — (r — i acos(#)) 4 ip4 . 

This formulation has several advantages: i) It is a first 
order gauge invariant description, ii) It does not rely 
on any frequency or multipole decomposition, iii) It can 
be used to evolve 3+1 dimensional spacetimes without 
any assumption about symmetries (to deal with the fi- 
nal stage of orbiting binary black holes), iv) The Weyl 
scalars are objects defined in the full nonlinear theory and 
it can be argued that evolving them with the linear theory 
provide a reliable description of the perturbations [25]. In 
addition, the Newman-Penrose formulation constitutes a 
simple and elegant framework to organize higher order 
perturbation [24]. 

The numerical integration of the linear Teukolsky 
equation in the time domain using Boyer-Lindquist co- 
ordinates is done closely following reference [26] . We use 



This allows us to reduce the dimensionality of the Teukol- 
sky equation from 3+1 to 2+1. In addition this decompo- 
sition into modes can be applied to the output of the full 
numerical code with the advantage of handling 2D fields 
instead of 3D ones. Typical evolutions of the Teukol- 
sky equation used a grid size of ng x n r » = 40 x 1200, 
with -18 < r*/M < 78 for signals of t - 100M, and we 
filled in initially with zeroes (or used extrapolations) the 
grid-points outside the full numerical domain. Finally, 
the computation of the energy and momenta radiated is 
performed using the formulae of Ref. [24], Sec. III.C. 

It worth stressing here that the Teukolsky equation can 
be written in any coordinate system. We are using it in 
the Boyer-Lindquist coordinates for present convenience, 
but if full numerical codes including excision black hole 
interiors turn out to be more practical in Kerr-Schild-likc 
coordinates, it may be convenient to evolve perturbations 
in a Kerr-Schild slices of the Kerr metric [27]. 



VII. APPLICATION TO A SINGLE ROTATING 
KERR HOLE 

We have done extensive testing of our method on var- 
ious toy models. In an earlier incarnation, we tried our 
approach successfully on axisymmetric head-on collision. 
As we have described, we have done a lot of work gener- 
alizing our method to include orbital cases with angular 
momentum on the final black hole. We have checked 
our equations explicitly on exact Boyer-Lindquist Kerr 
data, but in our real numerical simulations we will not 
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reproduce the Boyer-Lindquist coordinates exactly and 
is it useful to get some measure of how important the 
coordinate differences are for the radiation. Addition- 
ally, our calculation requires the use of four computer 
codes, a code for the numerical simulation (Cactus [28]), 
a specialized code (called Zorro) running within the sim- 
ulation which calculates the all the quantities needed for 
producing Cauchy data (see Appendix), a code (called 
TeukCauchy) which runs after the simulation defining 
the background black hole and constructing the Cauchy 
data from output of the evolution on various time slices, 
and the Tcukolsky evolution code (TeuiCode). The Kerr 
case has been a key source of rigorous tests of all these 
codes. 

We performed a highly nontrivial test of our set up by 
applying the entire procedure to Kerr initial data, evolv- 
ing a single Kerr hole full numerically, finding a "back- 
ground" black hole in the numerical data and defining 
a tetrad, extracting the Cauchy data and continuing its 
evolution with the Teukolsky equation. Ideally, the fi- 
nal result should produce no radiation. In practice, the 
computed radiation energy and waveforms will give us a 
measure of the error with which we can determine such 
quantities. 

This is a nontrivial test because the full numerical evo- 
lution is performed with vanishing shift and the singu- 
larity avoiding maximal slicing. This is in contrast to 
the Boyer-Lindquist lapse and the non-vanishing Boyer- 
Lindquist shift for rotation parameter a/M = 0.8 in the 
example shown in Fig. 8. In addition the Cauchy data 
for ^4 an d dtipi are computed with a tetrad adapted 
to the numerical spacetime, which must then be trans- 
formed according to Sec. V A in order to nearly reproduce 
the Kinnersley tetrad in the perturbative limit. These 
complications mean, in particular, that the data passed 
between our first three codes is not expected to be ap- 
proximately vanishing, but must sum to zero in the end. 
In practice our result is subject to both numerical error 
and false radiation caused by an inexact identification of 
the background coordinates and tetrad, perhaps produc- 
ing nontrivial Cauchy data which we then evolve via the 
Teukolsky equation. The results of the whole procedure 
are summarized in Fig. 8. The levels of spurious radi- 
ation are around 10 _5 M. Results after relatively short 
evolution times converge quadratically toward zero with 
increasing resolution. The longer evolutions are affected 
by the location of a close exterior boundary, and are im- 
proved when we move the boundary outwards by 50%. 
(As discussed above we will use much more distant outer 
boundaries for our astrophysical applications.) All this 
indicates that the coordinate effects are, so far, smaller 
than the numerical effects, which in turn tend to produce 
radiation about two orders of magnitude smaller than 
the radiation we are interested in. Notably these results 
are achieved with lower resolutions and closer boundaries 
than the typical resolutions of runs we performed for two 
black holes starting from the ISCO configuration used in 
Ref. [6]. 



*w(10 M), Kerr, a/M=0.8 



64x128* ,dx=M/8 
96x192? ,dx=M/8 
1 28x256?, dx=M/1 6 




T/M 

FIG. 8: The total radiated energy for evolved Kerr hole for 
different resolutions and boundary location. 



For the sake of completeness we mention two further 
tests which we performed for two black hole initial data: 
i) The mass scaling of the whole procedure. Since Ein- 
stein equations scale with the total ADM mass, we made 
a full numerical run with initial mass equal 2 and com- 
pare the scaling of the Cauchy data, post-processing and 
final waveforms with the mass equal 1 case. This proved 
to be a very useful test for the corresponding set of four 
codes we used to compute each of the above stages, ii) 
Reducing the initial separation of the holes from that of 
the ISCO to one quarter of it we reach the close limit 
regime and can compare with the results with the known 
analytic expressions [29] and scaling with the separation 
as well as angular dependence of ip^ and dtip4- 



VIII. DISCUSSION 

The Lazarus approach to binary black holes combines 
three treatments, each adapted to one of three stages of 
the dynamics, the far- limit, non- linear-interaction, and 
close-limit (one black hole) regimes. In this paper we 
have provided a detailed explanation of how numerical 
simulation and Teukolsky equation perturbative dynam- 
ics can be interfaced to provide a complete description 
of gravitational radiation arising from the post-orbital 
binary black hole dynamics. 

This technology makes it possible, for the first time, 
to apply numerical relativity, to the non-linear dynami- 
cal interaction of these systems. In our approach to this 
unknown regime, we have identified several parameter se- 
quences which make a connection to better-studied cases. 
An "L-sequence" , allows us to increase the separation 
from close-limit regime to ISCO, an "a-sequence" which 
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allows us to connect boosted head-on collisions studied 
in 2D to the ISCO with a fixed magnitude of each black 
hole's momentum, and a "P-sequence" through which wc 
connect to head-on collisions of resting black holes by 
varying the magnitude of the momentum to the ISCO 
value [6]. These important studies will give us some 
understanding of how the dynamics from ISCO config- 
urations relate to and differ from the simpler problems 
treated so far by numerical simulations and close-limit 
studies. 

After establishing such a basis for understanding the 
near ISCO regime dynamics we can reach out further, 
and seek to firmly establish the relation of these ISCO 
black hole configuration to astrophysics. As we have 
discussed in detail for the numerical simulation/close- 
limit interface, the dual approach to dynamics in overlap- 
ping validity regions provides a vital consistency check 
on the reliability of the results. A key goal, which we 
can now begin to approach, is to provide the same sort 
of consistency studies to the far-limit/numerical simu- 
lation interface and to thereby establish a firm astro- 
physical foundation for expensive, and difficult, numeri- 
cal work. Again, a handful of initial data sequences are 
appropriate for beginning to evaluate the connection of 
numerical/close-limit results to astrophysical problems. 
Within the effective-potential method which we have 
taken advantage of in order to define of ISCO data are a 
natural "Pi-sequence" of pre-ISCO stable circular orbits. 
Moving up this sequence toward more separated black 
holes asymptotically eliminates the features of these data 
which may be less astrophysical. Similarly, it is pos- 
sible to define curves through the parameter space of 
our initial data family which approaches the trajecto- 
ries defined by the Buonanno-Damour extension of the 
post-Newtonian method. The application of more ad- 
vanced numerical techniques [30] , which we are presently 
undertaking, should make it possible to begin generat- 
ing waveforms from farther up these sequences. Still, 
though, the effective potential method initial data se- 
quence is an imperfect stand-in for robust interface with 
the post-Newtonian method, which we expect to be ul- 
timately required. Since a primary concern about the 
astrophysical relevance of numerical/close- limit results is 
artificial radiation content in the initial data; another 
useful line of research is comparison studies of waveforms 
alternative initial data sets which would be equivalent in 
their astrophysical interpretation. These will provide a 
measure of significance of this interpretive indeterminacy 
to the gravitational radiation. Promising work with evo- 
lutions from Kerr-Schild-like initial data, for which an 
alternative instantiation of the Lazarus approach is un- 
der development [31], should enable an example of such 
comparative work. 

Another area of study which can now be pursued is to 
develop some preliminary indications of the effect spin 
has on the waveforms generated in the post-inspiral dy- 
namics. The effective potential approach provides a de- 
scription [32] of the effect small amounts of individual 



black hole spin have of the ISCO initial data. We are 
applying our approach at first instance to cases of spin 
parallel and antiparallel to the orbital angular momen- 
tum. 

Eventually numerical simulations will run routinely for 
hundreds or thousands of M, having begun from estab- 
lished astrophysical data. But the possibility for observa- 
tion is beginning almost immediately, and until now we 
have not met the needs of observers who express that any 
additional information about the the final stage of binary 
black holes may be extremely important [33] . Our efforts 
have shown that a crucial requirement for producing re- 
sults relevant to observers is to adapt numerical evolu- 
tions to astrophysical problems. Numerical relativity is 
ready, now, to begin answering questions about binary 
black holes in the near ISCO and pre-ISCO regime. 
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APPENDIX A: A PRACTICAL CONSTRUCTION 
OF THE CAUCHY DATA 

Here we describe a procedure for calculating tp4 which 
allows us to cleanly separate the specification of the back- 
ground black hole from the numerical simulation, as is 
practical for studying variations of the background met- 
ric and background coordinates. Since we have not yet 
determined the background black hole at the time of evo- 
lution we must compute a larger set of quantities which 
can then be transformed to the desired result after the 
background is specified. The steps are: 

(a) Compute the numerical 3 + 1 tetrad components as 
in Eq. (5.6) 

(b) With this tetrad, using the Cauchy data on a nu- 
merical time slice directly compute all five Weyl scalars 
ipo . . .ip4 and their time variations dgipo • ■ ■ do^A as fol- 
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lows: 



slice: 



ip4 — ~Rijkin l m : >n k m l 
+2Rojki(n°m3n k rh l — m°n J 'n k m l ) 
+RQjOi(n rh j n m l + m°n j m°n l - 2n°fh j fh°n l ) 

(Al) 

V>3 = Rijkil t n j fh k n l 

+R 0jkl (l°n^m k n l - n°m j l k n l - n°Pm k n l + fh°nn k n l ) 
+Ro j0 l{l n j m n l - l°n j n°m l - n°Pfh°n l + n°Pn°m l ) 

(A2) 

1P2 = Rijkil l m J rh k n l 

+R 0jkl (l°m j m k n l - n°rh j l k m l - m°l :i rh k n l + m°n j l k m l ) 
+Ro j0 i(l m j m n l - fm 3 n m l - m l j m°m l + n Pm°m l ) 

(A3) 

i>i = R ij kin i l j m k l l 

+Ro jk i(n°l j m k l l ~ l°m^n k l l - l°n 1 m k l l + m°Pn k l l ) 
+R 0j oi(n Pm°l l - n°Pl°m l - ZV'm¥ + l°nH a m l ) 

(A4) 

ipa = Rijkil l m?l k m l 
+2Ro jkl (l°m j l k m l - m°l j l k m l ) 
+Ro j0 i(l m j l m l + m°l j m°l l - 2l°m j m°l l ) 

(A5) 



do ^4 = doRi j kin l m 3 n k fh l 
+2d R 0j ki(n°rh j n k m l - m°n j n k m l ) 
+d R 0j oi{n fh 3 n°fh l + fh°n 3 m°n l - 2n°m : >fh n l ) 

(A7) 

do^3 = d R i jkil t n : >rh k n l 

+d Roj ki{l°n j rh k n l - n°rh j l k n l - n°Prh k n l + m°n j l k n l ) 
+d Ro j oi(l n i fh°n l - l°n j n°m l - n°Pm°n l + n°l j n°m l ) 

(A8) 

dofo = d R ij kil t m : >fh k n l 

+d R ojk i(l°m j m k n l - n°m j l k m l - m°l j m k n l + fh°nH k m l ) 
+doRr m (l m j m°n l - l°m j n°m l ~ m° P fh° fh l + n°l j m°fh l ) 

(A9) 

do^i = d R ij kin l l i m k l l 

+d R 0j ki{n°Pm k l l - l°m^n k l l - l°n j m k l l + m°Pn k l l ) 
+d R 0j oi(n l i m l l - n°l j l°m l - l°n j m°l l + l°nn°m l ) 

(A10) 

dotpo = d Ri j kil t m 3 l k m l 

+2d Q RQ 3 ki{l°mn k m l - m°l j l k m l ) 

+doRojoi (l°m j l°m l + m°l j m°l l - 2l°m j m°l l ). 

(All) 



The derivatives involved in the above expressions can be 
where R a bcd is the four-dimensional Riemann tensor: computed in terms of the data on the Cauchy hypcrsur- 

facc using Einstein's equations, 



Rijkl = Rijkl + 2Ki[ k K[]j 



Rfli 



Ojkl 



K 



j[k,i] 



Vi ]P 



Rojoi — Rji — KjpKf + KKj 



(A6) 



We compute the time variations only from data on the 



doRijkl = -47V {Ki[ k Ri]j - Kj[ k Ri]i 

~2 R ( K ^[k9i]j - Kj[k9l]i) } 

+^9i[kdoRt\ 3 - 2g j [ k d a Ri] i - gi[ k gi]jd R 
+2Ki[ k doKijj — 2Kj[ k doKqi, 



d Rojki = doKj^q + d T p j[k K l]p + T p j[k d K l]p , 



(A12) 

p> 
(A13) 



doRojoi — 



d Rii - 2Kf,d K 



~2NK JP KPK? + KjtdoK + Kd Q Kji\ (A14) 
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where 

d K = NK pq K pq - V 2 iV, 

(A15) 

daKij = N [Rij + KK l3 - 2K ip K p j - N^ViVjN] . 

(A16) 

d R = 2KP c <d Q K pq + AN K vq K v a K sq - 2K d a K 

(A17) 
(A18) 

dolt- = -2V {l (NK 3) k ) + V k {NK l0 ). (A19) 

These time variations may be precisely the background 
time derivatives we want if have already specified the 
background spacetime and can set N = AT( ), the back- 
ground lapse. On the other hand when we want to de- 
termine the background independently of the numerical 
simulation code we can produce the required informa- 
tion for later processing from all five quantities with the 
choice N = 1. 

(c) Next we define the background coordinate system 
as described in Section IV. This allows us now to refer to 
quantities defined in the background coordinates. If we 
have used N = 1 in constructing the time variations we 
must next translate the time variations calculated above 
to genuine background time derivatives with a set of cor- 
rections for the effect of the background lapse and shift. 
For the lapse, 

<9 Vo — » ad ipo + 2 a,f ipo + 2 a§ ipi (A20a) 
ad tp! +a ti .^i + ]-ag (ifo + 3?A 2 j(A20b) 



doi>i 
d tp3 



ad ip 2 + a e (ipi + ips) (A20c) 
ad ip 3 - a. f ^3 + ^ atf {ip 4 + 3^2lA20d) 
adoip4 — 2a.fip4 + 2a gip 3 (A20c) 



where atf and a g are related to the lapse for the back- 
ground Kerr metric, 



ot a — 



= i/^N(O),0- 



(A21) 



The shift corrections are 
Bo^Po -*d iPo +iP%A + iP%i>i+N k d k ^ (A22a) 



(V>2 - i>2) + N k d k ^ (A22b) 
5oV-2 +^^^i+^)+ Nk dki>2 (A22c) 

4v-3 +\ p$. (v>4 + + (Vs - 

(V>2 - + N k d k yj 3 (A22d) 
d V>4 ^d 4>4 ~iP%i>4+iP%^+N k d k ^4 (A22e) 



£>o4>2 
do*p3 



where and ff- are related to the shift for the back- 
ground Kerr metric, 



g rr JV (0),r 



1% = 



(0),fl 



(A23) 



(0),' 



(d) 



The Weyl scalars corresponding to the transformed 
tetrad defined in Eq. (5.9) , ip4 and d t tp4 can then 
be respectively expressed as a linear combination of 
the five numerical-tetrad Weyl scalars tp 4 . . . ipo and the 
Boip4 ■ ■ ■ doipo, as given in Eqs.(A22), with coefficients de- 
pending on the background coordinates r and 8, and on 
M and a: 



1p4 



AF\F% 



[(Va* + 1 - i)Vo 



+MA(y/A? + l - l)V>i - 6A 2 i/j 2 

-4iA(VA 2 + 1 + 1)V> 3 + (VA 2 + 1 + 1) 2 V 4 ] (A24) 



d t ip4 = 



AF\Fl 



[(VA 2 + 1 - l) 2 9 ^o 



+AiA{^A 2 + 1 - l)5bVi - 6A 2 4^2 
-4L4(v / I 2 Tl + 1)5 ^3 + (v 7 ^ 2 + 1 + l) 2 5oV^4 



(e) 



We use e lm ^ KS decomposition which is affected by the 
ip transformation given in Eqs. (4.4) and (6.4). This 
transformation is implemented at the end of the calcula- 
tion by ip4 — > e lm¥>o//sct ?/>4 and likewise for d t ip4- Note 
that the calculation of the last shift correction term in 
Eqs. (A22) can be also conveniently carried over after 
step (e) rather than in step (d). 
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